A Mortar BDD method for solving flow in stochastic discrete fracture networks

In geological media, the large variety and complex configurations of fractured networks make it difficult to describe them precisely. A relevant approach is to model them as Discrete Fracture Networks (DFN)[10, 19], with statistical properties in agreement with in situ experiments [15, 13, 14]. A DFN is a 3D domain made of 2D fractures intersecting each other. Steady state flow in DFN is considered, the rock matrix is assumed impervious. Following a Monte-Carlo approach, a large number of DFN has to be generated and for each, a flow problem has to be solved whatever the complexity of the generated networks. Moreover time and memory costs for each simulation should be as lower as possible. A nonconforming discretization of DFN allows to reduce the number of unknowns and facilitate mesh refinement. Sharp angles are managed by a staircaselike discretizations of the fractures’ contours [34]. The non-matching feature at the fractures’ intersections is handled via a Mortar method [4, 5, 1] developed for DFN in [33, 34] for a mixed hybrid finite element formulation. It consists in defining, for each intersection between fractures, master and slave sides. Due to the staircaselike discretizations, a shared edge may be labeled several times with master and/or slave properties, it is called in the paper a multi-labeled edge. Continuity conditions are enforced between the unknowns on both sides. The derived linear system has only inner and master traces of hydraulic head as unknowns. The matrix A of this system is a symmetric definite positive (SPD) arrow matrix in presence of Dirichlet boundary conditions [34]. The challenge is to solve such linear systems with millions of unknowns [17]. Direct solvers (like Cholmod [11]) are very efficient for small systems but suffer from a high need of RAM memory when the system size becomes too large. Among iterative solvers, multigrid methods are very efficient for most networks but for some, the convergence rate is very slow [35, 17]. Preconditioned Conjugate Gradient (PCG) is efficient and robust for every network tested [35]. The natural decomposition of the matrix A in subdomains encourages the use of domain decomposition methods [7, 36, 31, 24]. The Schur complement of the matrix A is SPD and yields an interface system with only master unknowns. This interface system can be solved iteratively